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Abstract 

Spike-timing dependent plasticity (STDP) is an organizing principle of biological neural net- 
works. While synchronous firing of neurons is considered to be an important functional block 
in the brain, how STDP shapes neural networks possibly toward synchrony is not entirely clear. 
We examine relations between STDP and synchronous firing in spontaneously firing neural pop- 
ulations. Using coupled heterogeneous phase oscillators placed on initial networks, we show nu- 
merically that STDP prunes some synapses and promotes formation of a feedforward network. 
Eventually a pacemaker, which is the neuron with the fastest inherent frequency in our numerical 
simulations, emerges at the root of the feedforward network. In each oscillatory cycle, a packet 
of neural activity is propagated from the pacemaker to downstream neurons along layers of the 
feedforward network. This event occurs above a clear-cut threshold value of the initial synaptic 
weight. Below the threshold, neurons are self-organized into separate clusters each of which is a 
feedforward network. 
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I. INTRODUCTION 



Synchronous firing of neurons lias been widely observed and is considered to be a neural 
code that adds to firing rates. For example, experimental evidence suggests the relevance 
of synchronous firing in stimulus encoding feature binding , and selective attention 
[s], Q]. Collective dynamical states of neurons including synchrony may appear as a result of 
self-organization based on synaptic plasticity. Modification of synaptic weights (i.e., weights 
of edges in the network terminology) often occurs in a manner sensitive to relative spike 
timing of presynaptic and postsynaptic neurons, which is called spike-timing-dependent 
plasticity (STDP). In the commonly found asymmetric STDP, which we consider in this 
work, long-term potentiation (LTP) occurs when presynaptic firing precedes postsynaptic 
firingby tens of milliseconds or less, and long-term depression (LTD) occurs in the opposite 
case p]. The amount of plasticity is larger when the difference in the presynaptic spike time 
and the postsynaptic spike time is smaller [5|. 

The asymmetric STDP reinforces causal pairs of presynaptic and postsynaptic spikes and 
eliminate other pairs. Based on this property of STDP, how STDP may lead to various forms 
of synchronous firing has been studied in both experiments and theory. Synchronous firing 
in the sense of simultaneity of spike timing can be established in recurrent neural networks 
when the strength of LTP and that of LTD are nearly balanced jol . Large-scale numerical 
simulations suggest that reproducible spatiotemporal patterns of spike trains self-organize in 
heterogeneous recurrent neural networks [3, Q] • Self-organization of clusters of synchronously 
firing neurons that excite each other in a cyclic manner has also been reported j^, [l^ . 

We previously showed that STDP leads to formation of feedforward networks and entrain- 
ment when there is a pacemaker in the initial network [11]. We considered random networks 
of coupled oscillators whose synaptic weights change slowly via STDP. We assumed that 
the oscillators have a common inherent frequency except a single pacemaker whose inherent 
frequency is larger. By definition, the rhythm of the pacemaker is not affected by those of 
other oscillators. The network generated via STDP is a feedforward network whose root 
is the pacemaker. In a final network, a spike packet travels from the pacemaker to the 
other neurons in a laminar manner. The neurons directly postsynaptic to the pacemaker 
fire more or less synchronously just after the pacemaker does. These neurons form the first 
layer. These neurons induce synchronous firing of the neurons directly postsynaptic to them, 
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which define the second layer. In this fashion, a spike packet starting from the pacemaker 
reaches the most downstream neurons within relatively short time, which resembles the phe- 



nomenology of the synfire chain [1^. Compared to the case of frozen synaptic weights, a 
pacemaker entrains the rest of the network more easily with STDP in the meaning that 
entrainment occurs with smaller initial synaptic weights. 

The previous work does not explain how pacemakers emerge. No matter whether the 
pacemakers are intrinsic oscillators or network oscillators, they pace rhythms of other ele- 
ments without being crucially affected by other rhythms. Although some pacemakers may 
be "robust" oscillators whose rhythms are insensitive to general input, a more natural ex- 
planation may be that pacemakers emerge through synaptic plasticity in a neural network 
in which pacemakers are initially absent. In this case, emergent pacemakers do not have to 
be robust oscillators; their rhythms can change in response to external input. The emergent 
network topology makes such neurons pacemakers by eliminating incoming synapses. A 
neuron would fire with its own rhythm if it is not downstream to any neuron. This scenario 

nn 

is actually the case for two-neuron networks [HI, 113| . Here we are concerned to networks of 
more than two neurons. An associated question is which oscillator may become a pacemaker. 

In this work, we numerically investigate recurrent networks of coupled phase oscillators 
subject to STDP. We show that, when the initial synaptic weights are strong enough, STDP 
indeed yields feedforward networks so that downstream neurons are entrained by an emergent 
pacemaker. To our numerical evidence, the emergent pacemaker is always the neuron with 
the largest intrinsic frequency. Below the threshold for entrainment, STDP leads to the 
segregation of the initial neural network into subnetworks of feedforward networks. 



II. MODEL 



A. Coupled phase oscillators 

We model dynamics of neural networks by coupled phase oscillators whose synaptic 
weights are plastic. Although a majority of real neurons fire in the excitable (i.e., fiuctuation- 
driven) regime, for tractability we use phase oscillators, which fire in an oscillatory manner. 
Generally speaking, phase transitions are more easily and clearly determined in the oscilla- 
tory regime than in the excitable regime. This is a reason why collective neural dynamics 
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Q,Q including ones associated with STDP flfl Il6| have been analyzed in the oscillatory 
regime actually to give insights into dynamics of neural networks possibly operating in the 
excitable regime. In the following, we report numerical results for = 3 and = 100. 

The state of neuron i {1 < i < N) is represented by a phase variable (pi G [0,27r). 
We identify 0j = and 0j = 27r. When 0j crosses in the positive direction, neuron i 
is defined to fire. We denote by tj and ti the spike time of presynaptic and postsynaptic 
neurons. If (pi crosses in the positive direction as time advances from t to t + At, we set 
ti = t + [2tt — (pi{t)]/[27i + (pi{t + At) — (pi{t)]At. As the initial condition, we set 0, = 
(1 < i < A^) for A^ = 3. We adopt this artificial initial condition to draw phase diagrams 
to systematically understand possible routes to synchrony via STDP. For A^ = 100, (pi{0) is 
picked randomly and independently for each i from the uniform density on [0, 27r). Neuron i 
is endowed with inherent frequency uJi so that it fires regularly at rate uJi/27r when isolated. 
Connectivity between neurons is unidirectional and weighted, consistent with the properties 
of chemical synapses. The set of edges in a network is denoted by E. In other words, 
(j, i) E E if neuron j is presynaptic to neuron i. Dynamics of the coupled phase oscillators 
are given by 

^ = + 4t- Y1 3x1 sin(0j - (pr) + aii, (1) 

where {k) is the average in-degree of neuron z, gji is a synaptic weight, and E,i represents the 
standard Gaussian white noise independent for different i. As a result of the phase reduction 
theory 171], the coupling term in the oscillatory regime is generally given by a 27r-per iodic 
function of the phase difference (pj — (pi under the assumption of weak coupling. This is also 
the case for pulse coupling, for which averaging an original pulse coupling term over one 



oscillatory cycle results in a coupling term as a function of (pj — (pi [ij]. Modeling realistic 
synaptic coupling needs a coupling term that contains higher harmonics [14]. However, our 
objective in the present paper is not to precisely describe the neural dynamics but to clarify 
general consequences of STDP under the oscillatory condition. We thus employ the simplest 
coupling term (i.e sinusoidal coupling). 

For A^ = 3, we set the amplitude of the noise a = 0.0071 so that the phase transitions 
are sharp enough and artificial resonance that is prone to occur when inherent frequencies 
satisfy MiUi = Mjujj for small integers Mj and Mj (1 < ^ < j < A^) is avoided. Accordingly, 
an independent normal variable with mean and standard deviation a^/At = 0.00071 is 
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added to each neuron every time step At = 0.01; we use the Euler-Maruyama integration 
scheme with unit time At. To determine the phase transitions for = 100, we do not 
apply dynamical noise because, up to our numerical efforts, the numerical results do not 
significantly suffer from artificial resonance. In some other simulations with = 100, we 
add different amplitudes of dynamical noise to examine the robustness of the results. 



B. STDP 



With STDP, gji is repeatedly updated depending on spike timing of neuron j and i. 
Specifically, LTP occurs when a postsynaptic neuron fires slightly after a presynaptic neuron 

n 

does, and LTD occurs in the opposite case [5|. We assume that synaptic plasticity operates 
much more slowly than firing dynamics. We denote by and A~ the maximum amount of 
LTP and that of LTD incurred by a single STDP event. Most of previous theoretical work 



supposes that A is somewhat, but not too much^ 
rates and to keep neurons firing [gI, Q, y, [ol, lol . 



arger than A'^, to avoid explosion in firing 



Ul]. Therefore we set = 0.9. How a 

single spike pair specifically modifies the synaptic weight is under investigation [s], [isl , and 
triplets or higher-order combination of presynaptic and postsynaptic spikes rather than a 
single presynaptic and postsynaptic spike pair may induce STDP [19]. However, we consider 
the simplest situation in which STDP modifies synaptic weights in an additive manner and 
the amount of STDP is determined by the relative timing of a presynaptic and postsynaptic 
spike pair. A single synaptic modification Agji triggered by a spike pair is represented by 

v4+exp(-^) ti-U<0 

where r is the characteristic time scale of the learning window, which is known in experiments 
to be 10-20 ms sj. Given that inherent frequencies of many pyramidal neurons roughly range 
between 5 and 20 Hz, r is several times smaller than a characteristic average interspike 



interval. Therefore, following 



llj, we set r = 1/6 x 27r/ci;, where u; is a typical value of spike 



frequency that is used to determine Ui. Following our previous work [ll|], we set oo = 8.1. 
Because learning is slow compared to neural dynamics, A~ must be by far smaller than a 
typical value of g. To satisfy this condition, we set A~ = 0.001 for A^ = 3. When A^ = 100, 
average in-degree (k) is set equal to 10. This implies that a neuron receives about five to 
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ten times more synapses than when = 3. To normahze this factor, we set it A = 0.0001 
for = 100. 

We assume that gji is confined in [0, gmax]] all the synapses are assumed to be excitatory, 
because the asymmetric STDP explained in Sec. [T] has been found mostly in excitatory 
synapses. Because dynamical noise is assumed not to be large, all the synaptic weights 
usually develop until gji almost reaches either Qmax or 0, until when we run each simulation 
run. Note that, even if gji = is reached, still belongs to E. The upper limit gmax is 
determined so that a notion of synchronization that we define in Sec. Ill CI does not occur 
when gji = gmax, V(j,i) G E. Accordingly, we set g^ax = 7.5 and g^ax = 15 for A^ = 3 and 
A^ = 100, respectively. 



C. Measurement of synchrony 

To obtain the threshold for synchrony in Sees. IIII Al and IITIB 11 we start numerical simula- 
tions with the initial condition gji = gQ, V(j, i) G E. There are various notions of synchrony. 
We focus on the possibility of frequency synchrony in which neurons fire at the same rate. 
In the oscillatory regime, frequency synchrony is commonly achieved in two main ways. One 
is when neurons are connected by sufficiently strong mutual coupling. Then they oscillate at 
the same rate and with proximate phases. The other is when some neurons entrain others. 
When upstream neurons, which serve as pacemakers, entrain downstream neurons so that 
they are synchronized in frequency, synchronous firing in the sense of spike timing may be 



le same level in the hierarchy 
1201 ]. We explore possible 
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missing due to synaptic delay. However, neurons located at t 
relative to the pacemakers tend to have close spike timing 
emergence of such dynamics when pacemakers are initially absent in networks. 
We quantify the degree of frequency synchrony by order parameter r defined by 

2 



log 
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(3) 



where Ui = d(j)i/dt is the actual instantaneous frequency of neuron i when coupled to other 
neurons. If all the neurons fire exactly at the same rate, r would become negative infinity. 
In the actual frequency synchrony, r takes a large negative value mainly because of time 
discretization. We manually set = —4 for A^ = 3 and Vc = —9 for A^ = 100, so that r < 
corresponds to the full frequency synchrony. The value of fc for A^ = 100 is smaller than 



for = 3 for two reasons. First, in the numerical simulations determining the degree of 
frequency synchrony, dynamical noise is present for = 3 and absent for = 100. Second, 
we are concerned to the frequency synchrony of a// the neurons so that V^Sj' '^j')^ 

is small regardless of A^; we have to normalize the prefactor 1/A^ in Eq. 



III. RESULTS 



A. Networks of three neurons 



Our goal is to understand dynamics of large neural networks. As a starting point, we ex- 
amine network evolution and possibility of frequency synchrony using small networks, which 
will help us understand dynamics of large networks. Two-neuron networks were previously 



analyzed 



We need at least three neurons to understand competition between differ- 



ent synapses, pruning of synapses, and effects of heterogeneity. Accordingly, we examine 
dynamics of different three-neuron networks under STDP. 



1. Complete graph 

Consider the complete graph [Fig. [T](a)], in which every pair of neurons is bidirectionally 
connected. The complete graph does not survive STDP because LTP of a synapse implies 
LTD of the synapse in the reversed direction and the amount of LTD is assumed to be 
larger than that of LTP for the same time lag. We examine which synapses survive and 
whether frequency synchrony emerges through STDP. If a predetermined pacemaker exists 
in a network, the activity of the other neurons will be entrained into the rhythm of the 
pacemaker with sufficiently large initial synaptic weights, which was previously shown for 



N = 2 and A^ = 100 llj. Here we consider A^ = 3 and compare numerical results when 
a pacemaker is initially present and absent in the complete graph. Note that the effective 
initial network when the pacemaker neuron 1 is initially present is the one shown in Fig.[T](b), 
because the synapses toward the pacemaker are defined to be entirely ineffective. 

First, we examine the relation between heterogeneity in inherent frequencies, initial 
synaptic weights, and synchrony. We expect that small heterogeneity and large initial 
synaptic weights favor synchrony. To focus on phase transitions, we reduce the number 
of parameters by setting all the initial synaptic weights equal to Qq and restrain inherent 



frequencies ui, 0^2, and (cui > > ^^3) by imposing, ui — UJ2 = UJ2 — UJ3 = Au, where 
UJ2 = 8.1. Numerically obtained phase diagrams are shown in Figs. [2t^a) and [2](b) for the 
cases in which a pacemaker is initially present and absent, respectively. The results are 
qualitatively the same for the two situations. The neurons get disconnected and fire inde- 
pendently as a result of STDP for sufficiently small qq or sufficiently large Au (blue regions 
labeled D). A feedforward network whose root is the fastest oscillator emerges for sufficiently 
large qq or sufficiently small Auj (yellow regions labeled A). Then all the neurons rotate at 
frequency Ui. In the intermediate regime (green regions labeled C), final synaptic weights 
satisfy (723 ~ dmax and gi2, gn, g2i, gsi, and (732 ~ 0. In this case, neuron 2 entrains neuron 
3 so that they oscillate at frequency uj2, whereas neuron 1 gets disconnected and oscillates 
at frequency Ui. We rarely observed the case in which neuron 1 entrains 2 (or 3) and neuron 
3 (or 2) gets isolated. Although ui — uj2 = uj2 — UJ3, neuron 1 is more likely to segregate 
from the network than neuron 3 is. Quantitatively speaking. Figs. [2](a) and [2](b) indicate 
that the entrainment of the entire network by the fastest neuron (i.e., neuron 1) is to some 
extent easier to realize when the pacemaker is initially absent than present (yellow regions 
labeled A). In Figs. [2](a) and[2]^b), the phase diagrams are disturbed along vertical lines at 
Auj ^ 2.7. This artifact comes from the fact that Ui, uj2, and 003 approximately satisfy the 
resonance condition (i.e., MiUi = M2UJ2 = M3UJ3 with small integers Mi, M2, and M3). In 
some of the following figures, similar disturbance appears along special lines. We can wash 
away these artifacts by increasing the amount of dynamical noise. However, we prefer not 
doing so to prevent the boundaries between different phases from being blurred too much. 

Next, to examine what happens when ui, 002, and uj^ change independently, we set go = 
0.15, UJ2 = 8.1, and vary Aui = uoi — UJ2 and Auj2 = UJ2 — uo^. Numerical results with and 
without a pacemaker are shown in Figs. Ws) and[2](d), respectively. Figures [2](c) and[2]^d) 
are similar to each other, except yellow spots in the red region (labeled B) in Fig.[2]^c). These 
spots represent entrainment facilitated due to the artificial resonance condition satisfied by 
uji,uj2i and ^^3. Both in Figs. Ws) and[2|^d), (723 is easier to survive than gi2 is, consistent with 
Figs. [2](a) and[2]^b). This is indicated by the fact that the phase of the frequency synchrony 
of the three neurons (yellow regions labeled A) extends to a larger value of Auj2 > along 
the line Aui = than to the value of Aui > along the line Auj2 = 0, and that the phase 
in which neuron 2 entrains 3 (green, C) survives up to a larger value of Auj2 than the value 
of AuJi up to which neuron 1 entrains neuron 2 but not neuron 3 (red, B). 
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To examine the cause of the asymmetry in Figs. [2]^c) and[2]^d) along the two hnes Aui = 
and Auj2 = 0, we analyze a two- neuron network with asymmetric initial synaptic weights 
shown in Fig. [2]^a). The two neurons h and / have inherent frequency Uh and Ui (< Uh). The 
weights of the synapse from neuron h to neuron / and that from neuron / to neuron h are 
denoted hj Qf and gt, respectively. When Aui = and Auj2 > in the three- neuron network, 
neurons 1 and 2 are synchronized almost from the beginning, in both frequency and phase, 
because ui = uj2- This is true if a trivial condition gi2+g2i > is satisfied. Then the network 
is reduced to the two- neuron network by identifying ojh = oji = uj2, oji = u^, g f = gi^, + g23, 
and gb = {g2i + fi'3i)/2. When, Aui > and Auj2 = in the three-neuron network, neurons 
2 and 3 are synchronized in frequency and phase as far as (723 + 932 > 0. Then the network 
is reduced to the two-neuron network with ujh = uji, ui = uj2 = uj^, gf = {gi2 + gi3)/2, and 
gb = g2i + gsi- For these two situations, we calculate the threshold for frequency synchrony 



in the two-neuron network using the semi-analytical method developed in Because all 
the synaptic weights are initially equal to go in Fig. [2], the initial condition for the two- 
neuron network is {gf,gb) = (2fi'o,5'o) for Aui = 0, Auj2 = Au > 0, and {gf,gb) = (fi'o,25^o) 
for Auji = Auo > 0, Auj2 = 0. The phase-transition curves for the frequency synchrony are 
shown in Fig. [3t^b), indicating that the threshold is larger along the Auj2 = line than along 
the AuJi = line. This is consistent with the three-neuron results shown in Figs. [2]^c) and 

Hd). 



2. Feedforward loop 

Other three-neuron networks, particularly feedforward ones, are presumably embedded 
in larger neural networks in the course of network evolution. First, we consider the network 
shown Fig. Il](a) as the initial network. 

Figure 111(a) is the phase diagram in which we vary Au = Ui — 002 = 002 — UJ3 and 5^0 = 
9i2 = 913 = 923- The original network shown in Fig. |4](a) survives STDP when initial 
synaptic weights are large or the heterogeneity is small (yellow region labeled A). In the 
opposite situation, all the neurons get disconnected and fire independently (blue, D). Neuron 
1 detaches from the network and neuron 2 entrains neuron 3 in the intermediate regime 
(green, C). 

The phase diagram in the Aco'i-Au;2 parameter space with go = 0.15 is shown in Fig.|l](b), 
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which looks similar to Figs. [2](c) and[2](d). As in the case of the complete graph, the situation 
in which neuron 1 entrains neuron 2 with neuron 3 isolated is less likely to arise than that 
in which neuron 2 entrains neuron 3 with neuron 1 isolated. 

3. Fan-in network 

Next, we examine dynamics starting from the fan-in network shown in Fig. [5]^a). In this 
network, neuron 3 is postsynaptic to two pacemaker neurons 1 and 2. We are concerned to 
which neuron entrains neuron 3. 

First, we examine the case in which two synapses are initially equally strong and the 
inherent frequencies of the two upstream neurons are different. Accordingly we set gi^ = 
923 = do, uji — Us = Aui, UJ2 — uj^ = Auj2, Qo = 0.2, and = 8.1. Figures [5](b) and 
Et^c) are the phase diagrams in the Au;i-Aco'2 space, with Fig. [5]^c) being an enlargement 
of Fig. [5]^b). There are principally four phases: neither neuron 1 or 2 entrains neuron 3 
(blue regions labeled D), both neurons 1 and 2 entrain neuron 3 (yellow. A), only neuron 
1 entrains neuron 3 (red, B), and only neuron 2 entrains neuron 3 (green, C). The phase 
diagram is symmetric with respect to the diagonal line Aui = Auj2. When uji and uj2 are 
too far from uo^, all the neurons get disconnected (blue, D). Both gi^ and (y'23 survive only 
when uji ~ 002 (yellow. A). This phase extends to the disconnection phase (blue, D) on the 
diagonal because, on this line, the firing of neuron 1 elicits LTP of both synapses so does 
firing of neuron 2. However, this situation is not generic in that uoi and UJ2 must be very close 
for this to happen. When uji and UJ2 are not close to each other and not too far from uj^, 
which upstream neuron entrains neuron 3 is not obvious. Figure [5](b) tells that a necessary 
condition for an upstream neuron to entrain neuron 3 is that the difference between its 
inherent frequency and c^s is less than ^ 1.0. This condition roughly corresponds to the 
requirement for the entrainment in the two-neuron feedforward network with = 0.2. This 
explains the two rectangular regions Acji > 1.0, Auj2 < 1.0, and AuJi < 1.0, Auj2 > 1.0 of 
Fig. [5]^b). In the remaining region (i.e. Ac^i < 1.0 and Auj2 < 1-0), the upstream neuron 
whose inherent frequency is closer to uj^, equivalently, the slower upstream neuron, largely 
wins the competition (regions marked by □). The faster upstream neuron entrains neuron 
3 when the inherent frequency of the slower upstream neuron is very close to (regions 
marked by Q)- The total size of the latter regions is much smaller than that of the former 
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regions. 

Starting with asymmetric synaptic weights, that is, gis 7^ (723, the upstream neuron more 
strongly connected to neuron 3 may entrain neuron 3. To investigate the interplay of this 
effect and heterogeneity in the inherent frequency, we perform another set of numerical 
simulations with ui = uj^ + 1, UJ2 = uJi + Au, gis = go, and (?23 = 90 + Ag^. The asymmetry 
in the initial synaptic weight is parameterized by AgQ. Figures [5](d)-[5]^f) shows the phase 
diagrams in the Au-Ago space for three different values of ui. On the singular line Auj = 
(i.e., uji = UJ2), Ago > 0, both upstream neurons entrain neuron 3. On the line Auj > (i.e., 
uji < UJ2), Ago = 0, neuron 1, whose inherent frequency Ui is closer to UJ3 than 002 is, entrains 
neuron 3 if Ui is not too apart from [Fig. E](d)]. This is consistent with the results in 
Figs. [S](b) andO^c). However, if g23 is sufficiently larger than gi^, neuron 2 overcomes the 
disadvantageous situation 002 — 003 > ui — to win against neuron 1 and entrains neuron 3. 
We confirmed that neuron 2 exclusively entrains neuron 3 when Au < and Ago > (not 
shown) . 

B. Networks of many neurons 

In this section, we use networks of heterogeneous = 100 neurons to examine what net- 
work structure and dynamics self-organize via STDP when we start from random neural net- 
works. The inherent frequencies of the neurons are independently picked from the truncated 
Gaussian distribution with mean 8.1, standard deviation 0.5, and support Ui G [7.6,8.6]. 
We assume that every neuron has (k) = 10 randomly selected presynaptic neurons on av- 
erage so that an arbitrary pair of neurons is connected by a directed edge with probability 
{k)/{N — 1) ^ 0.1. Except in Sec. IIIIB 3[ where we investigate effects of heterogeneity, the 
initial synaptic weight is assumed to be go common for all the synapses. We vary go as a 
control parameter. 

1. Threshold for frequency synchrony 

We compare how STDP affects the possibility of entrainment and formation of feed- 
forward networks when a pacemaker is present and when absent. To this end, we fix a 
random network and a realization of tUj (1 < < A^). Without loss of generality, we assume 
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TABLE I: Comparison of the threshold for frequency synchrony and the actual mean frequency 
of the neurons (lD) in the frequency synchrony. We calculated (lD) by averaging the instantaneous 
frequency over all the neurons and over the last ten unit times of the simulation. 





Pacemaker 


Present 


Absent 


STDP 


Present 


gc = 9.8 

(cD) = 8.60 


gc = 0.72 
(lD) = 8.60 


Absent 


9c = 51 
(w) = 8.60 


gc = 0.93 
(a)) = 8.08 



uji > UJ2 > ■ ■ ■ > ojn- For the network with a pacemaker, we make the fastest neuron a 
pacemaker. By definition, the rhythm of the pacemaker is not affected by those of the other 
neurons even though the pacemaker is postsynaptic to approximately {k) neurons. Using 
the bisection method, we determine the threshold value of above which all the neurons 
will synchronize in frequency. 

The results without dynamical noise (i.e., cr = 0) are summarized in Table I. When 
the pacemaker is present from the beginning, STDP drastically reduces the threshold for 



entrainment After entrainment, all the neurons rotate at the inherent frequency of 

the pacemaker, that is, uoi = 8.60. When a pacemaker is initially absent, STDP reduces 
the threshold for frequency synchrony by 34%. Facilitation of frequency synchrony in the 
absence of the initial pacemaker is consistent with the results for the complete graph with 
= 3 (Fig. [2]). In this situation, the scenario to frequency synchrony is different between 
the presence and the absence of STDP. With STDP, the fastest oscillator eventually entrains 
the entire network when the initial synaptic weight is above the threshold, as in the case of 
the network with a prescribed pacemaker. Without STDP, the fastest oscillator does not 
entrain the other neurons. The realized mean frequency 8.08 is close to the mean inherent 
frequency of the 100 neurons. This suggests that frequency synchrony in this case is achieved 
by mutual interaction, rather than by one-way interaction underlying the entrainment by the 
fastest neuron. Therefore, in networks without predetermined pacemakers, STDP enables 
emergence of pacemakers and changes the collective dynamics drastically. 
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2. Network dynamics 



For 0" = 0, example rastergrams when there is initially no pacemaker and go = 1.0, which 
is above the threshold value 0.72 (see Table I), are shown in Fig. O Figures [Hl^a) andlHI^b) 
correspond to the initial and final stages of a simulation run under STDP, respectively; 
frequency synchrony appears as a result of STDP. Figure [6]^c), which is an enlargement of 
Fig. [6](b), shows that the fastest neuron entrains the other neurons and that faster neurons 
tend to fire earlier in a cycle. Figure [7] shows the time course of the degree of synchrony 
r. Around t = 1.2 x 10^, r sharply drops, and all the neurons start to oscillate at the 
same frequency. The effective network defined by the surviving synapses in the final state is 
drawn in Fig. [HI The neurons are placed so that the horizontal position represents relative 
spike time in a cycle. With this ordering, the neurons form a feedforward network. In other 
words, after STDP, if a presynaptic neuron fires later than a postsynaptic neuron in a cycle, 
this synapse is not present. 

Partial entrainment occurs when is slightly or moderately smaller than the threshold 
value 0.72. Circles and crosses in Fig. [3 represent the actual frequency after transient and 
the inherent frequency of the each neuron, respectively, when g^ = 0.5. The neurons with 
the same actual frequency belong to the same cluster. Each cluster forms a feedforward 
network emanating from an emergent pacemaker. Figure [9] indicates that the neurons are 
divided into two clusters and one isolated neuron. Neuron 2 entrains 85 other neurons all 
of which are slower than neuron 2, neuron 6 entrains 12 slower neurons, and neuron 1 gets 
isolated. In this and further numerical simulations we performed, the root of a feedforward 
subnetwork is always occupied by the fastest neuron in the cluster. 

Whether two neurons eventually belong to the same cluster is determined by where these 
neurons are located on the initial random network and by how close their inherent frequencies 
are. If go is smaller than the value used for Fig. [9l two neurons have to be closer in Ui to stay 
connected after STDP. Then the number of clusters increases, and the number of neurons 
in a cluster decreases on average. 
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3. Robustness against dynamical noise and heterogeneity 



To examine the robustness of the numerical results reported in Sec. HUB 2t we perform 
additional numerical simulations with dynamical noise and random initial synaptic weights. 
We draw initial gji {j, i) & E randomly and independently from the uniform density on 
[0,25(0], where go = 1.0. 

With a = 0.081, the rastergram and the actual frequency of the neurons after transient 
are shown in Figs. [TOT a) andflOlfb). respectively. With a = 0.081, the standard deviation of 
the accumulated noise in a unit time, which is equal to cr, corresponds to 1% of the phase 
advancement estimated by the mean inherent frequency of the oscillators, which is equal to 
8.1. The rastergram [Fig. [TUlfa)] is indicative of full entrainment. Indeed, all the neurons 
eventually rotate at the inherent frequency of the fastest neuron [Fig. fTOT b)]. With a = 0.405, 
the neurons are divided into six synchronous clusters of size 31, 26, 19, 12, 5, 4, plus 
three isolated neurons [Figs. [TOT c) and [TOT d)]. With cr = 0.81, many neurons, particularly 
faster ones, rotate at their inherent frequencies [Figs. [TOTe) andfTOVf)]. Consequently, there 
are many clusters of neurons. The frequency synchrony within each cluster is blurred by 
dynamical noise. 

In sum, emergence of entrainment via STDP survives some dynamical noise and het- 
erogeneity in the initial synaptic weights. We have confirmed that, when the entrainment 
occurs, it is quickly established at around t = 10^, and the fastest oscillator is located at 
the root of the feedforward network, as in Fig. [8l 



4- Network motifs 

We investigated the evolution of three-neuron networks in Sec. IIII Al because we expect 
that these results have something common with evolution of such subnetworks in large 
networks. The results in Sec. IIII Al predict the following: 

• Bidirectional edges do not survive STDP, and feedforward networks of size three will 
be relatively abundant after STDP. Subnetworks abundant in a large network relative 
to the case of random networks with the same mean degree (or other order parameters) 



are called network motifs 



22']. The hypothesis that feedforward networks are motifs in 



large neural networks is consistent with the observations in C. elegans neural networks 
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22|. 

• As a result of STDP, a neuron has at most one effective upstream neuron unless 
multiple upstream neurons are very close in frequency. 

There are 13 connected network patterns of three nodes. How often each pattern appears 
in a network with = 100, relative to the random network, can be quantified by the 
Z score [22]. The Z score is the normalized number of a pattern in the network, where 
normalization is given by the mean and the standard deviation of the count of the pattern 
based on independent samples of the randomized networks. A pattern with a large Z score 
is a motif of the network with A^ = 100. 

Figure [H] shows the Z score of each pattern before (circle) and after (square) STDP, 



calculated by m-finder 23j. We set a = (i.e., no dynamical noise) in this analysis. The error 
bar shows a range of one standard deviation based on ten simulation runs in each of which 
we draw a different initial random network and a different realization of cuj (1 < z < 100). 
Before STDP, the neural network is a directed random graph, so that the Z score for each 
pattern is around zero, meaning that no pattern is overrepresented or underrepresented 
significantly. After STDP, the feedforward network whose emergence and survival were 
observed in Sec. IIII Al (i.e., pattern 5 in Fig. [TT!) and patterns consistent with this (i.e., 
patterns 1 and 2) are overrepresented. These are motifs of our final networks. Pattern 4 
is also a motif in spite of our negation in Sec. IIII Al because the two upstream neurons in 
pattern 4 have the same actual frequency. They are generally different in inherent frequency 
but share a more upstream ancestor. As the example network in Fig. [H] shows, existence of 
multiple paths from a neuron to another due to branching and uniting of edges is compatible 
with STDP. The other network patterns are not significant or underrepresented. These 
results are further evidence that feedforward networks are formed by STDP in heterogeneous 
neural networks. 



IV. DISCUSSION 

We have shown using heterogeneous coupled phase oscillators that feedforward networks 
spontaneously emerge via STDP when the initial synaptic weights are above the threshold 
value. When this is the case, the pacemaker emerges at the root of the feedforward network 
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and entrains the others to oscillate at the inherent frequency of the pacemaker. Although 



these results have been known for two-neuron networks [11|, [ISj , we have shown them for the 
cases of three and more neurons and quantified the phase transitions separating frequency 
synchrony and asynchrony. The route to frequency synchrony is distinct from a conventional 
route to frequency synchrony that occurs when mutual, but not one-way, coupling between 
oscillators is strong enough. Some results obtained in this work are unique to the networks 
without a prescribed pacemaker. First, the emergent pacemaker is the fastest oscillator 
neuron according to our extensive numerical simulations. Note that all the oscillators fire 
at this frequency in the entrained state, whereas they fire at the mean inherent frequency 
of the oscillators when the frequency synchrony is realized by strong mutual coupling in the 
absence of STDP. Second, when the initial coupling strength is subthreshold, the neurons 
are segregated into clusters of feedforward networks. Third, our numerical evidence suggests 
that entrainment under STDP occurs more easily when a prescribed pacemaker is absent 
than present. 

In spite of a wealth of evidence that real neural circuits are full of recurrent connectiv- 
ity 24]) feedforward structure may be embedded in recurrent neural networks for reliable 



transmission of information 
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25j. Feedforward transmission of synchronous volleys in 



rather homogeneous neural networks as those used in this work serves as a basis of repro- 
ducible transmission of more complex spatiotemporal spike patterns in more heterogeneous 
networks. Such patterns may code external input or appear as a result of neural compu- 
tation 0, Feedforward structure is also a viable mechanism for traveling waves often 
found in the brain 261]. Although computational roles of feedforward network structure are 
not sufficiently identified, our results give a support to the biological relevance of feedfor- 
ward networks. The formation of feedforward networks, which we have shown for oscillatory 
neurons, is consistent with numerical results for more realistic excitable neurons subject to 
STDP 271]. The neurons that directly receive external input may be more excited and fire 
at a higher rate compared to other parts of a neural circuit. Our results suggest that such 
a neuron or an ensemble of neurons is capable of recruiting other neurons into entrainment 
and creating feedforward structure. 

We assumed the additive STDP with the nearest-neighbor rule in which the dependence of 
the amount of plasticity on the current synaptic weight and the effects of distant presynaptic 
and postsynaptic spike pairs, triplets, and so on, are neglected. Generally speaking, evolution 
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of synaptic weights is affected by the implementation of the STDP rule [3]. However, we 
believe that our results are robust in the variation of the STDP rule as far as it respects the 
enhancement of causal relationships between presynaptic and postsynaptic pairs of neurons. 
Our preliminary numerical data with excitable neuron models suggest that the results are 
similar between the multiplicative rule |l8| and the additive rule [Kato and Ikeguchi (private 
communication)]. Recent reports claim the relevance of acausal spike pairs in the presence 



of synaptic delay 
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281]. This and other factors, such as different time scales of LTP and 



LTD j27j, may let bidirectional synapses survive as observed in in vitro experiments |29 |. 
Incorporating these factors is an important future problem. 

We have ignored inhibitory neurons for two reasons. First, our main goal is to identify 
phase transitions regarding synchrony with a simple model. Second, specific rules of STDP 
are not established for inhibitory neurons, albeit some pioneering results 30|. Taking in- 



hibitory neurons into account, preferably in the subthreshold regime, warrants for future 
work. 
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FIG. 1: Complete graph (a) without a pacemaker and (b) with a pacemaker. 
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FIG. 2: (Color) Phase diagrams for the complete graph in the [(a) and (b)] Auj-go and [(c) and 
(d)] AuJi-Auj2 spaces. One pacemaker neuron is initially present [(a) and (c)] or absent [(b) and 
(d)]. We run numerical simulations 20 times for each pair of parameter values. We add the red 
element of the RGB color scheme by the maximum amount divided by 20 when gi2 survives in 
a simulation run. Similarly, the green is added when §23 survives, and the blue is added when 
all the neurons get disconnected. Yellow regions appear when both gi2 and 523 survive, since the 
combination of red and green is yellow. In this case, it turns out that gis also survives. We verified 
that no other connectivity, such as survival of gi^ without survival of (712 or g23, and survival of g2i, 
g^i, or (^32, appears except at points near phase transitions and resonance. Near phase transitions, 
we exclude such exceptional runs from the statistics. In the resonance regions (e.g., Au ~ 2.7 and 
go ~ 0.4), the three neurons may remain connected. In this situation, however, synaptic weights 
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FIG. 3: (a) Two-neuron network, (b) Threshold for frequency synchrony for the two-neuron 
networks corresponding to the Atui = hne and the Aijj2 = hne in Figs. W^c) and[2jd). 




FIG. 4: (Color) (a) Feedforward loop, [(b) and (c)] Phase diagrams for the feedforward loop in 
two different parameter spaces. See Fig. [2] for the color code. 



(b) 



(c) 




(f) 




(e) 



0.0 I 

0.0 0.2 0.4 0.6 




12 3 




12 3 

Ago 








B 




C 


■ 















FIG. 5: (Color) (a) Fan-in network, [(b) and (c)] Phase diagrams for the fan-in network in the 
Aa;i-Aa;2 space, with (c) being an enlargement of (b). We set go = 0.2. [(d), (e), and (f)] Phase 
diagrams in the Auj-Ago space. We set = 0-2, (d) ui = uj^ + 0.8, (e) uji = uj^ + 1.0, and (f) 

LVi = UJ3 + 1.2. 




FIG. 6: Rastergrams for (a) initial 5 time units and (b) final 5 time units of a simulation run. (c) 
is an enlargement of (b). We set N = 100, go = 1.0, and a = 0. The neurons are aligned according 
to the order of the inherent frequency. 




FIG. 7: Time course of the degree of synchrony when N = 100, go = 1.0, and a = 0. The values 
of r are plotted every 10000 time units. 




FIG. 8: Final network structure when 
Pajek 3l- 



= 100, go = 1.0, and a = 0. The network is drawn by 
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FIG. 9: Segregation into clusters when N = 100, = 0-5, and a = 0. Inherent frequencies (+) 
and actual frequencies after STDP (o) are shown. We estimate the actual frequencies from the 
phase shifts with bins of width 10 time units. 
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FIG. 10: Results for N = 100, go = 1.0, with dynamical noise and heterogeneity in the initial 
synaptic weights. We set [(a) and (b)] a = 0.081, [(c) and (d)] a = 0.405, and [(e) and (f)] 
a = 0.81. Rastergrams for 5 time units after transient are shown in (a), (c), and (e). The neurons 
are aligned according to the order of the inherent frequency. Inherent frequencies (+) and actual 
frequencies (o) are shown in (b), (d), and (f). Because of the dynamical noise, we estimate the 
actual frequencies from the phase shifts with bins of width 10^ time units. 
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FIG. 11: Normalized abundance of different three-neuron network patterns. We set = 100, 
go = 5.0, and a = 0. Circles and squares correspond to the initial and final stages of the simulation 
runs, respectively. 



